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^ i Abstract 

• I describe a method, particularly suitable to implementation by com- 

Q^ . puter algebra, for the derivation of low-dimensional models of dynamical 

I systems. The method is systematic and is based upon centre manifold 

^~>' theory. Computer code for the algorithm is relatively simple, robust and 
flexible. The method is applied to two examples: one a straightforward 



O I pitchfork bifurcation, and one being the dynamics of thin fluid films. 



1 Introduction 

Deterministic evolution equations, either ordinary differential equations or partial 
differential equations, are used to describe dynamics in the physical world. Often, 
these equations allow many modes of behaviour that are of little physical interest 
in particular applications, for example sound waves are neglected in many appli- 
cations of fluid mechanics. The essential dynamical behaviour of the system is 
then determined by the evolution of a subset of the possible modes, the "critical" 
modes. Rapid oscillations or heavy damping characterise the modes that need 
to be eliminated from consideration. If the amplitudes of the modes are viewed 
as co-ordinates in a state space then the dynamical behaviour of the system cor- 
responds to motion along some trajectory. When many of the modes are heav- 
ily damped, trajectories are rapidly attracted to some low- dimensional invariant 
manifold, which may be parameterised by the amplitudes of the critical modes. 
This geometric picture is the heart of the application of centre manifolds Q 
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to the rational construction of low-dimensional model systems by the elimina- 
tion of physically uninteresting fast modes of behaviour. Such low-dimensional 
dynamical models are significantly easier to analyse, simulate and understand. 

Applications of the techniques have ranged over, for example, triple convection 
[0, feedback control 0, economic theory 0, shear dispersion [18|, nonlinear 



oscillations [^, beam theory [Q, flow reactors 0, and the dynamics of thin fluid 



films . New insights given by the centre manifold picture enable one to not 



only derive the dynamical models, but also to provide accurate initial conditions 



[pil , [T^, boundary conditions [^, and the treatment of forcing p. 

In the application of centre manifold theory, the typical approach is to express 
the model explicitly in terms of asymptotic sums [§, |l], |2^, H, |4|, ^ e.g.]. 



To reduce the dynamics onto the centre manifold, one then has to substitute 
the asymptotic sums into the governing equations, reorder the summations, rear- 
range to extract dominant terms, and evaluate the expressions. While perfectly 
acceptable when done correctly, it leads to formidable working which obscures 
the construction of a model. Further, such asymptotic expansions, in common 
with the method of multiple scales, reinforce the notion that careful balancing of 
the "order" of small effects are necessary in the construction of a model rather 
than in its use in some situation. Instead, I propose in §2 an iterative method, 
based upon the residuals of the governing differential equations, for the construc- 
tion of such low-dimensional, dynamical models. The evaluation of the residuals 
is a routine algebraic task which may be easily done using computer algebra, as 
seen in the examples of §§2 and 3, by simply coding the governing differential 
equations; it replaces the whole messy detail of the manipulation of asymptotic 
expansions ( e.g. §5.4]). The aim of this proposed approach is to minimise 
human time by using a novel algorithm which can be simply and reliably im- 
plemented in computer algebra with relatively small inefficiencies in the use of 
computer resources. After all: "It is unworthy of excellent persons to lose hours 
like slaves in the labour of calculation" . . . Gottfried Wilhelm von Leibniz. 

The use of centre manifold theory in the construction of low-dimensional 
modes relies on three theorems. These specifically address dynamical systems 
written in the form: 

X = Ax + f(x,y), 
y = 5y + g(x,y), 

where: the overdot denotes d/dt; x(t) is m-dimensional; y{t) is n-dimensional 
(more generally, of infinite dimension); the eigenvalues of A have zero real-part; 
the eigenvalues of B have strictly negative real-parts bounded away from 0, < 
—7 < 0; and the nonlinear functions f and g are smooth and are at least quadratic 
near the origin. Then p4-5], §2] or [0, p5&p35]: 



Existence There exists a smooth m-dimensional centre manifold for ([l|) of the 
form y = h(x), tangent to y = at the origin. The dynamics on the centre 
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manifold are governed by 

a = Aa + f(a,h(a)). (2) 

That is, provided it has a correct spectrum of eigenvalues, a smooth dynam- 
ical system possesses a centre manifold with low- dimensional, self-contained 
dynamics. 

Relevance Let (x(t),y(t)) be a solution of the parent system (|l]) with initial 
point (x(0),y(0)) sufficiently small (in practice, sufficiently small may be 
quite wide), then provided the zero solution is stable, there exists a solution 
of (0) a(t) such that as t ^ cxd, 

x(t) = a(t) + O (e-^*) , y{t) = h(a(t)) + O (e-^*) , (3) 

where 7 > is some constant. 

That is, from a wide range of initial conditions, all solutions tend exponen- 
tially quickly to a solution on the centre manifold, and hence the dynamical 
system on the m- dimensional centre manifold, faithfully models the 
original. 

Approximation Seeking an approximation to the centre manifold, y = h(x), 
just substitute into the governing equations ([|), and use the chain rule to 
deduce that, given 



Mh = — 
ax 



Ax + f (x, h(x))] -Bh-g (x, h(x)) , (4) 



we wish to solve A^h = 0. Suppose that as x — >• 0, A^h = C(|x|''), 
then h(x) = h(x) + Odxl"^). It is often convenient to appeal to a more 
general assertion that explicitly accounts for constant parameters, say e: 
if Mh = O (|x|«, |en, then h(x, e) = h(x, e) + O (|x|^ jeH. (This more 
general form is particularly relevant to unfolding bifurcations and to the 



long- wave, slowly- varying approximation ||20|| .) 

That is, provided we can satisfy the governing equations to some order of 
accuracy, then the centre manifold will have been found to the same order 
of accuracy. 

Simple applications of these theorems may be found in the book by CarrQ]. 
Here I show how to use an algorithm, eminently suitable for computer algebra 
and based upon these theorems, to derive low-dimensional models of dynamical 
systems. In §2, I develop the formalism in general and in a bifurcation example. 
In §3 I show how straightforward it is to apply the techniques to a much more 
difficult problem, namely the flow of a thin fllm of fluid upon a solid substrate. 
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One important feature of the analysis is that we deal here with the problems 
in terms of the physical differential equations as given, and not in the abstract 
form A linear change of basis, such as (x, y) = Pu for some linear transfor- 
mation P, is all that is needed to transform from a physical description, in terms 
of physical variables u, into a form for a direct application of theory. 



2 Modelling a pitchfork bifurcation 

Consider the following variation to Burger's equation featuring growth, (1 + e)u, 
nonlinearity, uu^, and dissipation, u^x'- 

+ + + „,(0) = .W = 0. (5) 

for some function u{x, t). View this as an infinite dimensional dynamical system, 
the state space being the set of all functions u{x) on [0, vr]. 

For all values of the parameter e there is a fixed point at the origin, that 
is, a trivial equilibrium state is m = 0. A linearisation of the equation about 
this equilibrium, namely Ut = {1 + e)u + Uxx, has modes sinkx with associated 
eigenvalues = 1 — fc^ + e for wavenumbers k = 1, 2, . . .. Thus the k = 1 mode, 
sinx, loses stability as e crosses zero, and the system undergoes a bifurcation. 

To find the details of this pitchfork bifurcation is a simple task for low- 
dimensional modelling. Linearly, exactly at critical, e = 0, all modes decay 
exponentially quickly except for the critical mode sin x; it has a zero decay rate 
and therefore is long lasting; by the first theorem we are assured that there ex- 
ists a centre manifold. Nonlinearly, and for e and u{x) near 0, all modes decay 
exponentially except for the critical modes which have a slow evolution. Thus, 
exponentially quickly we can model the dynamics solely in terms of the evolu- 
tion of the amplitude of the sinx mode; I define a to be its amplitude. By the 
relevance theorem, the evolution of a in time forms an accurate one- dimensional 
model of the original infinite-dimensional dynamical system (BI). 



2.1 The iteration scheme 

1 now proceed to simultaneously develop a novel and powerful algorithm to derive 
such a low-dimensional model while applying it to the specific example (|^). In 
general, I address dynamical systems in the form 

ii = £u + f (u, e) , (6) 

where: 

• u{t) is the evolving state "vector" (corresponding to (x, y) in (|l|)), u{x,t) 
in the example; 



2: Modelling a pitchfork bifurcation 



A J Roberts, 22 March, 1996 



5 



• £ is a linear operator whose spectrum, as required by centre manifold the- 
ory, is discrete and separates into eigenvalues of zero real-part, the critical 
eigenvalues (corresponding to the modes x in (0)), and eigenvalues with 
strictly negative real-part (corresponding to the modes y in (|^)); in the 
example Cu = u + Uxx (with implicit boundary conditions); 

• e is a parameter, potentially a vector of parameters; 

• and f is a function which is strictly nonlinear when considered as a function 
in u and e together, that is, f is quadratic or higher order in u and e as 
they tend to 0, f = ew + uu^ in the example. 

For simplicity, I only treat the case where the critical eigenvalues of C are exactly 
0; let the multiplicity be m. Cases where the critical eigenvalues have a non-zero 
imaginary part, as in a Hopf bifurcation, may be handled with the same ideas as 
described herein, but with more complicating detail. The aim is to find a low- 
dimensional model a = g(a), such as (0), for the evolution of m "amplitudes" 
a. These low- dimensional dynamics occur on the exponentially attractive centre 
manifold which may be described parametrically as u = v(a). 

The first stage is to identify the m critical modes, that is, those associated with 
the eigenvalue zero; these are necessary in order to project the linear dynamics 
and nonlinear perturbations onto the slow modes of interest. They may be found 
from the nontrivial solutions, Vj, of £v = 0; in general we need the critical 
eigenspace and so may need to find all the generalised eigen-modes. Then, in 
terms of modal amplitudes Oj, a linear approximation to the centre manifold and 
the evolution thereon is simply 

u(^) ~ '^^j'^j — such that a ^ Qa, (7) 
j 

where V = [v^] and where Q may be chosen in Jordan form in the case of gen- 
eralised eigenvectors (in (|lD the m critical modes are x, and the linear approxi- 
mation is u = (x, 0) such that x = Ax). In the example (|^), the eigenvalue zero 
is of multiplicity m = 1 and so there is only the one critical mode, v{x) = sinx; 
hence the linear approximation is 

u{x, t) ^ a sin x such that d ~ . (8) 

To model the nonlinear dynamics, this linear approximation needs to be mod- 
ified by nonlinear terms; such modification would be equivalent to seeking the 
nonlinear shape of the centre manifold, y = h(x), in (^. 

The second stage is to seek iterative improvements to a given level of descrip- 
tion of the centre manifold and the low- dimensional evolution thereon. The aim 
is to find a low-dimensional description which satisfies the nonlinear dynamical 
equation (|^). As in Newton's method for finding the zero of a function, we use 



2: Modelling a pitchfork bifurcation 



A J Roberts, 22 March, 1996 



6 



the residual of the governing equations in order to guide corrections. The iter- 
ation scheme is successful as long as it ultimately drives the residual to zero to 
the desired order of accuracy — see the approximation theorem. Suppose that at 
one stage of the iteration we have the approximate model 

u ^ v(a) such that a ^ g(a) ; 

approximate because the residual of the governing differential equation (|) 

ii-Cu- f (u, e) = ^g-Cif- f (v, e) = (a", e*") , (9) 

for some order of error, q and r, and where a denotes |a|. We seek to find "small" 
corrections, indicated by primes, so that 

u ^ v(a) + v'(a) such that a ^ g(a) + g'(a) , 

is a better approximation to the centre manifold and the evolution thereon. The 
aim of each iteration is to improve the order of the errors {q and r) so that, by 
the approximation theorem, we improve the accuracy of the model. Substituting 
into the governing differential equation (^, and using the chain rule for time 
derivatives, leads to 



da da 



g') = £v + CV + f (v 



Given that it is impossible to solve this for the perfect corrections in one step, 
seek an approximate equation for the corrections of O (a^, e"^) by: 

• ignoring products of corrections (primed quantities) because they will be 
small, O (a^^ + e^''), compared with the dominant effect of the linear cor- 
rection terms; 

• and replacing tilde quantities by their initial linear approximation wherever 
they are multiplied by a correction factor (introducing errors O (a"^^^, e''^^))— 
such approximation slows the iteration convergence to linear, as opposed 
to the quadratic convergence which would be otherwise obtained (if only it 
were practical). 

Thus we wish to solve 

<9v (9v' 

— g + Vg' + —ga = Cy + CV + f (v, e) . 

It is not obvious, but provided it is arranged so that Q is in Jordan form, as is 
often physically appealing, we may significantly simplify the algorithm by also 
neglecting the term at a cost of increasing the number of iterations needed 
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by a factor no more than m, the multiphcity of the zero eigenvalue of C. Thus, 

Oao ~ dt 



rearranging and recognising that = % by the chain rule, we solve 



/:v' - vg' = 1^ - /:v - f (v, 6) , (10) 

for the primed correction quantities. In the example, we seek to solve 

, d'^v' I . dv ,^ ^dv d'^v 

^ + IT^ - ^ smx = — - (1 + e)w - V- — , 

ox^ ot ox ox^ 

which in the first iteration from the linear approximation (P) is 

d^v' 1 

v' + — — g' sin x = —ea sin x sin 2x . 

ox^ 2 

The great advantage of this approach is that the right-hand side is simply 
the residual of the governing equation (^ evaluated at the current, tilde, ap- 
proximation; in essence, this residual is the quantity defined in (^). Thus at 
any iteration we just deal with physically meaningful expressions; all the compli- 
cated expansions and rearrangements of asymptotic expansions, as needed by the 



method of multiple scales (e.g. [15, §3.5]) or earlier methods to find the centre 
manifold (e.g. §5.4]), are absent. There is, of course, a cost and that lies in 
evaluating the residual (which potentially has enormous algebraic detail), much 
of which is repeated at every iteration. However, with the advent of computer 
algebra, all this detail may be left to the computer to perform — such mindless 
repetition is ideal for a computer — whereas all a human need concern themselves 
with is setting up the solution of (p!OD and not at all with the detailed algebraic 
machinations of asymptotic expansions. 

The main detail is to solve equations of the form 

i:v'-Vg' = r, (11) 

for some given residual r. Recognise that there are more unknowns than com- 
ponents in this equation; its solution is not unique. The freedom comes from 
the fact that we can parameterise the centre manifold via the amplitudes a in 
an almost arbitrary manner. The freedom can only be resolved by giving a pre- 
cise meaning to the m amplitudes a. Often one does define a to be precisely the 
modal amplitudes (as is done implicitly for (|I]) by seeking a centre manifold in the 
form (x, h(x))) in which case we seek corrections v' which are orthogonal to the 
generalised eigenvalues, Zj, of the adjoint of £. In the example, £ is self adjoint 
under the inner product (mi,M2) = f Jq UiU2dx, and so the adjoint eigenvector is 
also simply z{x) = sinx and so I require {sinx,v') = 0. More general definitions, 
such as an energy related amplitude, give rise to similar considerations to those 
that follow. There are two approaches to solving (|11]). 
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1. Numerically, it is generally easiest to adjoin the amplitude condition to the 
equation and solve 



' C -V " 




" v' " 




r 















where Z = [zj]. 

2. However, algebraically it is usually more convenient to adopt the follow- 
ing procedure which is also familiar as part of other asymptotic methods. 
Rewrite ([TT|) as CV = Vg' + r and recognise that C is singular due to the 
zero eigenvalue of multiplicity m. We choose the m components of g' to 
place the right-hand side in the range of C; this is achieved by taking the 
inner product of the equation with the adjoint eigenvalues zj (this corre- 
sponds to considering just the x modes in (|1|)) and thus giving the set of m 
equations 

{Z,V)g' = -{Z,v). 

In the example, g' = —^jQsmxr{x)dx which in the first iteration gives 
g' = ea. This equation is known as the solvability condition. 

Having put the right-hand side in the range of C we solve CV = r = 
Vg' -|- r for v', making the solution unique by accounting for the definition 
of the amplitudes a. In the example, we solve the boundary value problem 
v' + v'xx = 'f'ix) such that v'{0) = v'{tt) = and that v has no sin a; 
component. In the first iteration, the problem v' + v'xx = ~\(^^ sin 2a; with 
the above conditions has the solution v' = ^a? sin 2a;. 

D 

Then the last step of each iteration is to update the approximations for the 
centre manifold shape and the evolution thereon. For the example, after the 
first iteration we deduce u ~ a sin a; -|- |a^sin2x, which shows the nonlinear 
steepening/flattening of negative /positive slopes, and that the evolution is d ~ ea, 
which exhibits the loss of stability of the fixed point a = as e becomes positive. 
Further iterations in the example lead to the centre manifold being given by 

M = a sin a; H ^"^^ H ^^"^^ 3x + O {a^, e^) , (12) 

on which the system evolves according to 

1 — - 

a = ea -^a^ + O (a^ e^) . (13) 

The relevance theorem assures us that this 1-D model of the original infinite- 
dimensional dynamical system (|]), is valid exponentially quickly in time. From 
the model (p!3D, for example, we deduce the quantitative shape of the pitchfork 

bifurcation: there are stable fixed points at a ~ ^J7/(1^^7/3) . Physically, these 
fixed points represent a balance between the nonlinear steepening of the uUx 
term, and the dissipation of Uxx- 
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2.2 Computer algebra implementation 

A principal reason for adopting this approach is because it is simply and reliably 
implemented in computer algebra. Based upon the above derivation, the general 
outline of the algorithm is: 

1. preliminaries; 

2. initial linear approximation; 

3. repeat until residual is small enough; 

(a) compute residual, 

(b) find solvability condition, 

(c) compute correction to the centre manifold, 

(d) update approximations. 

Complete details of a reduce program for the particular example @ follows. 
The reason for using reduceQ is that it has excellent pattern matching and 
replacement capabilities through its operator and let statements. 

% Find pitchfork bifurcation in u_t=(l+eps)u+uu_x+u_xx 
y„ a(t) measures amplitude of sin(x) component in u(x,t) 

y. 

on div; off allfac; on revpri; % improves appearance of output 
let sin(~x)*cos(~y) => (sin(x+y)+sin(x-y) ) /2 ; % a trig rule 

y. 

% Define the inverse operator of u+u_xx 
operator linv; linear linv; 
let linv(sin(~k*x) ,x) => sin(k*x)/(l-k"2) ; 
y Define inner product with sin(x) 
operator sindot; linear sindot; 

let { sindot(sin(x) ,x) => 1, sindot (sin(~k*x) ,x) => }; 
% 

depend a,t; % asserts that a depends upon time t 

let df(a,t) => g; % so da/dt is replaced by current g(a,eps) 

y. 

u:=a*sin(x); g:=0; % initial approximation 

y 

y iterate until PDE is satisfied (to requisite order) 

let ■[eps~2=0, a"4=0}; y discard high-order terms in a & eps 

repeat begin 



^At the time of writing, information about reduce was available from Anthony C. Hearn, 
RAND, Santa Monica, CA 90407-2138, USA. E-mail: reduce@rand.org 
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write eqn:=df (u,t) -(l+eps)*u-u*df (u,x)-df (u,x,x) ; 

gd:=-sindot(eqn,x) ; 

write u:=u+linv(eqn+sin(x)*gd,x) ; 

write g:=g+gd; 
end until eqn=0; 
7. 

; end; 

Observe the how this implements the algorithm. 
1. The preliminaries do the following. 

• M improves the appearance of the printed output, whereas £5 tells 
REDUCE that we wish to linearise products of trigonometric functions. 

• £7-9 defines the operator linv to act as the inverse of C: 

— declaring it linear tells REDUCE to expand sums and products in 
the first argument and to only leave functions of the second argu- 
ment inside the operator, for example, linv(easina;+2a^sin2a;,x) 
is expanded to 

ealinv (sin x,x)+2a^ linv (sin 2a; ,x) ; 

— the let statement on £9 defines the action of the operator as 
the solution to v' + v'^^ = sin kx, namely v' = yrp- sin kx, the tilde 
before the k on the left-hand side matches any pattern (no action is 
defined for the singular case k = 1 because the pattern sin(~k*x) 
does not match sin(x), but any appearance of linv(sin(x) ,x) 
usefully signals an error). 

• £10-12 similarly defines sindot to be the inner product operator (sinx, •), 
the let statement now being a list, enclosed within braces, of evalua- 
tion rules. 

• £14-15 establishes that the variable a is to firstly depend upon time, as 
we use a as the time dependent amplitude in the model, and secondly 
that time derivatives of a, df (a,t), are to be replaced by the value of 
g, at the time of replacement, as g is to store the current approximate 
model evolution equation such as (O). 



2. £17 simply assigns the linear approximation (H) of the centre manifold to 
be the initial value of the variables u and g. 

3. £20-26 performs the iterations. 

• £20 controls the truncation of the asymptotic approximation. It gives 
a list of transformations which tell reduce to discard any factor in 
or higher and any factor in or higher; thus all expressions are 
computed to an error of O (e^, a^). 
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• £22 computes the residual; observe how it is a very direct translation 
of the governing equation (j^) into REDUCE symbols, a user of this 
approach only has to implement the governing equations, all the messy 
details of the asymptotic expansions are dealt with by the computer 
algebra engine. The iteration is repeated until this residual is zero, to 

• £23 is the solvability condition giving corrections to g. 

• £24-25 solves for the correction to u and updates the current approx- 
imation. 

The program could be used to derive high-order effects in a or e simply by chang- 
ing the discarding of factors on £20. Different systems with the same linear 
structure may be analysed simply by changing the computation of the residual 
on £22. 

By way of comparison. Rand & Armbruster's pp27-34] macsyma code 
for constructing the centre manifold of a finite-dimensional dynamical system, 
based upon power series manipulation, uses 53 lines of active code (although 
some are for input of the dynamical equations) whereas the above algorithm has 
only 16 lines of active code. 



3 Thin film fluid dynamics 

We now turn to the modelling of an important physical problem, that of the flow 
of a thin film of viscous fluid upon a solid substrate. Examples include the flow 
of rainwater on a road or windscreen or other draining problems 0, paint and 
coating flows [ER E^, and the flow of many protective biological fluids |]I2|. For 



simplicity we restrict attention to 2D fluid flow and seek a model for the evolution 
of the depth of the fluid film; because the film is thin expect that the vertical 
structure of the flow is relatively simple. I show how centre manifold theory can 
partially justify such a model, and how to derive the model using the algorithm 
developed in this paper. 

This fluid flow problem has many nonlinearities, as shown in Figure |I|: not 
only is the advection in the Navier-Stokes equation described by a nonlinear term, 
but also the thickness of the fluid film is to be found as part of the solution and 
its unknown location is another source of nonlinearity. When a fluid layer is thin 



then, as in dispersion in a pipe ||18[, the dominant dynamical processes occur 
along the thin film of fluid. Across the fluid film, viscosity acts quickly to damp 
almost all cross-film structure; the only long lasting dynamics are those of the 
relatively slow spread of the film along the substrate. In a nonlinear problem, this 
signature of interesting dynamics on a long time-scale along with uninteresting, 
quickly dissipated modes indicates that centre manifold theory may be used to 
create a low-dimensional model of the interesting dynamics. 
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atmospheric pressure 



y = vix,t) 



Navier-Stokes equation 

q(x,y,t) Pix,y,t) 



h 



solid, u = V = 



Figure 1: schematic diagram of a thin fluid film flowing down a solid bed. 



Assuming no longitudinal variations in x, a linear analysis of the equations 
shows that there is one critical mode in the cross-film dynamics associated with 
conservation of the fluid, all others decay due to viscosity. Consequently, it is 
natural to express the low-dimensional model in terms of the film thickness ri{x, t). 
Seeking solutions which vary slowly along the film, that is, assuming d/dx is 
a small parameter like e used previously, a centre manifold analysis creates an 
effective model of the dynamics. Through the algorithm described herein, we find 
the long-term evolution of long-wavelength modes is approximately described by 



drj 
dt 



ld_ 

3 dx 
d_ 

dx 



(14) 
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where dashes or roman numerals on r] denote d/ dx. The centre manifold theorems 
reasonably assure us that this model is indeed relevant to the long-term dynamics 
of thin films. However, this assurance is not yet rigorous because of deficiencies 
in the preconditions of current theorems.0 

The first line of the model (JT^) is the standard leading, or "lubrication," 
approximation to thin film dynamics; the second line may be viewed either as 
correction terms, or as terms indicative of the error in the leading approximation. 
An interesting diversion is attributable to the fact that although this is a nonlinear 
problem, conservation of fluid applies no matter how thick the fluid layer (for 
example, the right-hand side of (p!^ can always be written as a gradient). Thus 
the analysis is valid for arbitrarily large variations in the thickness of the film! 
just so long as the variations are sufficiently slow. Theorems on the existence 

^The major constraint upon a rigorous application of the theorems by Gallay pH is that 
we need to treat d/dx as a small parameter, whereas in principle it is unbounded. However, 
provided we confine attention to functions slowly-varying in x, then such derivatives would 
indeed be small. 
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and approximation of such global (in rj) centre manifolds are also given by Carr 
i PP31-32]. 

3.1 Governing equations 

As shown in Figure |I| we solve the Navier-Stokes and continuity equations within 
the Newtonian fluid for the velocity field q: 

dq ^ 1_ _2 

— + q-Vq = Vp + z/Vq, 

Ot p 

V-q = 0. 

For simplicity, restrict attention to two-dimensional flow taking place in the xy- 
plane: the x-axis is aligned along the solid bed of the flow; the |/-axis is per- 
pendicular. The viscous flow must stick to the solid bed to give the boundary 
condition 

q=0 ony = 0. (15) 

The surface of the fluid, y = rj{x,t), evolves with the flow. Because the free- 
surface is unknown we not only need two boundary conditions for the Navier- 
Stokes equations, we also need an extra boundary equation in order to be able 
to find 1]. The kinematic condition is that the fluid flow, as given by the velocity, 
q, must follow the free-surface as given hj y = ri{x,t): 

dri dr] 

- = .-n- ony = ,. (16) 

The other boundary conditions come from the forces acting across the free- 
surface. Above the fluid film we suppose there is a very light and essentially 
inviscid fluid such as air (one-thousandth the density of water). 

• Since air is inviscid, it cannot sustain any tangential stress across the sur- 
face, thus 

Since the density of air is very low, Bernoulli's equation asserts that fluid 
stress exerted normally across the surface has to be constant, say T„ = —pa, 
equal and opposite to air pressure which without loss of generality I take 
to be zero. However, the effect of surface tension is like that of an elastic 
membrane, it causes a pressure jump if the surface is curved: positive if the 
fluid surface is convex; negative if it is concave. For both the normal stress 
and surface tension to oppose atmospheric pressure 




(l + r/'2)p = 2/i 



dv ^ ,2du , fdu ^dv^ 
dy ^ dx ^ \dy dx 



a7]" 
V 1 + V 
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where a is the coefficient of surface tension. 

Now non-dimensionahse these equations, refering all scales to a as surface 
tension is taken to be the driving force. For a reference length, suppose that h 
is a characteristic thickness of the thin film as shown schematically in Figure |I[ 
Then non-dimensionalise by writing the equations with respect to: the reference 
length h; the reference time /ih/a; the reference velocity U = cr/ii, and the 
reference pressure a/h. With these choices, and in non-dimensional quantities, 
we solve the Navier-Stokes and continuity equations 



' dq 



V-q = 0, (19) 



where 7?. = ^ is effectively a Reynolds number characterising the importance 
of the inertial terms — it may be written as Uh/u for the above characteristic 
velocity. Observe that due to the absence of 71 in the model (p^, to its order 
of accuracy the inertial terms have no influence on the dynamics — it is creeping 
flow. The non-dimensional equations are subject to the bottom boundary condi- 
tion ([T5|) , and on the surface the kinematic condition ([I6D and the two dynamic 
conditions, (p!7| ) and 



(l + v") 



P 



dy ^ dx ^ 




7]" 

on y = rj. (20) 



TI 



In the Navier-Stokes equation I have neglected body forces. If the presence 
of gravity were to be acknowledged, then it would appear in the non-dimensional 
combination of the Bond number b = gph'^/a. This may be neglected if b is 
extremely small {h <^ 0.2cm for water). Including gravity, assuming b is small 
but not negligible, is not a signiflcant complication and I leave its inclusion as 
an exercise for the reader — perhaps the simplest interesting case is that of a fluid 
film on a vertical substrate. 



3.2 Linear picture 

The construction of a centre manifold model rests upon an understanding of the 
dynamics linearised about some fixed point (the spectrum of the linear operator 
is crucial in theory and application). Here the fixed point is a film of constant 
thickness, non-dimensionally 1, and of zero velocity and pressure. The dynamics 
linearised about this fixed point are based on velocities and rj — 1 being small. 
Thus, as well as neglecting products of small terms, boundary conditions on the 
free-surface are approximated by their evaluation on the approximate surface 
y = 1. The linearised equations are then 

= -Vp + v^q, 
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V-q = 0, 

q = on ?/ = 0, (21) 



V on y = 1, 



drj 
dt 

du dv ^ 

^ + ^ = on?/=l, 

oy ox 



dv 

p = 2— r]" on y = 1. 

dy 

To investigate the dynamics of "long waves" on this thin film we treat d/dx as 



small, that is slowly- varying in x, as justified formally in ||2^ and more rigorously 



for capillary-gravity waves in H. "Linearly" then we can base the construction of 



a centre manifold model on the limit when there are no longitudinal variations: 
^ = 0. Neglecting all x derivatives and seeking solutions proportional to c^* 
leads to the following. Firstly, Aq = is an eigenvalue associated with the mode 
7] = const, and u = v = p = 0. It is the presence of this eigenvalue which 
indicates the existence of a centre manifold model. Other modes exist, namely 
Un = sin7r(n— |)?/, Vn = Pn = Vn — ^ = with eigenvalues A„ = — 7r^(n — 1)^/7^ for 
n = 1,2, . . .. Thus for long- waves there exists one eigenvalue of the dynamics, 
whereas all the rest of the eigenvalues are strictly negative, bounded above by 
—7 ^ —2.5/7^ (the smaller the Reynolds number, the more creeping the flow, 
and the faster the decay of the non-critical modes). Hence, by the relevance 
theorem, expect the dynamics of thin films to exponentially quickly approach 
a low-dimensional centre manifold based on the mode corresponding to the 
eigenvalue, the film thickness 77, and characterised by slow variations in x. 



3.3 Iterative construction 

The first task is to decide how to parameterise the centre manifold model. Since 
the critical mode is rj = const., p = u = v = 0, it is appropriate to use the 
film thickness f] as the parameter. This is especially useful since rj has a direct 
physical meaning. Thus a linear description of the centre manifold is 

u = v = p = 0, (22) 

and 1] free to vary according to ^ ~ 0. 

Now we organise an iteration scheme to refine the description of the cen- 
tre manifold. The order of error will be characterised by the number of spatial 
derivatives of 77; for a slowly- varying function higher-order derivatives are asymp- 
totically smaller than lower-order derivatives. Suppose that the fields q^i], y) and 
p{ri, y) are an approximation to the centre manifold "shape," with the associ- 
ated approximate evolution ^ ~ .^('7)- Seek equations for corrections to this 
description by: 
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substituting into the governing equations:^ 

q = q(ry) + q'(r]) , 
P = P{r])+P{ri), 

such that — = g{r]) + g'{i]) ; 



dt 



• omit products of corrections; 

• omit x-derivatives of corrections as both corrections and are small; 

• approximate other terms involving corrections by replacing the current ap- 
proximation, tilde quantities, with the initial linear approximation (here 
zero); 

• rearrange the equations. 

The upshot is that we solve equations 



dp' 
dy'^ dy 
dv' 
dy 

with boundary conditions 



du 



n 



dt 

' dv 
V • q 



dx 

+ q ■ Vi] + ^-V^v 
J dy 



(23) 
(24) 
(25) 



p - 2 



du' 
dy 
dv' 
dy 



l + 7]'^]p + 2 



du dv 
dy dx 
dv 



J'2 



du 



dy dx 



' du ^ dv' 
dy dx 



rj" 



on y = 1] . 



q 

9 + 9' 



on y = , 

V + v' — (u + u')— on y = ri . 
dx 



(27) 
(28) 
(29) 



However, the unknown location of the free surface of the film is a major 
technical difficulty. One way to proceed is to scale the vertical coordinate, C = 
y/rj, so that the free surface corresponds to C = 1 precisely. Because rj varies 
with X and t, this scaling of y affects space-time derivatives and so plays havoc 



■^Be careful not to confuse corrections, indicated by primes, and derivatives of rj with respect 
to X, also indicated by primes. 
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with details of the governing equations. Under the change of coordinates from 

{x,y,t) to 

^ = x, c = y/vix^t), r = t, 

the chain rule shows that derivatives transform according to 

d d ^7]' d d d ^fj d d Id 
dx V 9C ' dt Or rj d( ' dy r] d( 

Fortunately we may implement the analysis in computer algebra and thus relegate 
virtually all such details to the computer. 

The only places where we need to explicitly consider these rules are in the 
terms in the equations for the corrections, u', v', and p'. However, these only 
involve y derivatives, see (^^- |27|) , which simply transform ~f^- Thus, we 

multiply the residuals on the right-hand sides of these equations by the appropri- 
ate power of rj, as seen in lines 40, 47 and 53 of the following reduce program. 

°/o Construct slowly-varying centre manifold of thin film fluids. 

°/o Allows for large changes in film thickness. 

% 

on div; off allfac; factor d,h,re; % improves printing 
°/„ use stretched coordinates: ys=y/h(x,t) , xs=x, ts=t 
depend xs,x,y,t; 
depend ys,x,y,t; 
depend ts,x,y,t; 

let i df(~a,x) => df (a,xs)-ys*h(l)/h(0)*df (a,ys) 

, df(~a,t) => df (a,ts)-ys*g/h(0)*df (a,ys) 

, df(~a,y) => df (a,ys)/h(0) 

, df(~a,x,2) => df (df (a,x) ,x) }; 
% solves -df (p,ys)=rhs s.t. sub(ys=l ,p)=0 
operator psolv; linear psolv; 
let ■[psolv(ys~~n,ys) => (1-ys" (n+l))/(n+l) 

,psolv(ys,ys) => (l-ys"2)/2 

,psolv(l,ys) => (1-ys) }; 
y„ solves df (u,ys,2)=rhs s.t. sub(ys=0,u)=0 & sub(ys=l ,df (u,y) )=0 
operator usolv; linear usolv; 

let {usolv(ys"~n,ys) => (ys~ (n+2)/(n+2)-ys)/(n+l) 

,usolv(ys ,ys) => (ys"3/3-ys)/2 

,usolv(l,ys) => (ys"2/2-ys) }; 
y„ use operator h(m) to denote df(h,x,m) 
operator h; 
depend h,xs,ts; 
let i df(h(~m),xs) => h(m+l) 

,df (h(~m) ,xs,2) => h(m+2) 
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,df (h(~ni) ,ts) => df(g,xs,m) }; 

y. 

% linear solution 
u:=0; v:=0; p:=0; g:=0; 
7. 

% Iteration. Use d to count the number of derivatives of x, 
% and throw away this order or higher in d/dx 
let d~7=0; 

curv:=h(2)*d"2*(l-h(l)~2*d"2/2+3*h(l)"4*d"4/8-5*h(l)"6*d~6/16) ; 
repeat begin 

7„ continuity & bed 

write ceq:=-df (u,x)*d-df (v,y) ; 

write v:=v+h(0)*int(ceq,ys) ; 

7„ vertical momentum & normal stress 

write veq:=re*( df (v,t)+u*df (v,x)*d+v*df (v,y) ) 

+df (p,y) -df (v,x,2)*d-2-df (v,y,2) ; 
write tn:= sub(ys=l , -p* (l+h(l) "2*d"2) +2*(df(v,y) 

+h(l)~2*d~3*df (u,x)-h(l)*d*(df (u,y)+df (v,x)*d)) 
-curv ) ; 

write p:=p+h(0)*psolv(veq,ys) +tn; 

% horizontal momentum & bed & tangential stress 

write ueq:=re*( df (u,t)+u*df (u,x)*d+v*df (u,y) ) 
+df (p , x) *d -df (u , X , 2) *d"2-df (u , y , 2) ; 

write tt:=-sub(ys=l, (l-d"2*h(l)"2)*(df (u,y)+df (v,x)*d) 
+2*h(l)*d*(df (v,y)-df (u,x)*d) ); 

write u:=u+h(0) "2*usolv(ueq,ys)+h(0)*tt*ys ; 

write g:=sub(ys=l,v-u*h(l)*d) ; 
end until (veq=0)and(tn=0)and(ueq=0)and(tt=0)and(ceq=0) ; 
; end; 

Recognise the algorithm in this code. 
1. Preliminaries 

• The change of coordinate rules are implemented in £5-12, where 
xs denotes ^, ys denotes (, and ts denotes r. Then wherever in the 
algebraic details we need these rules they are automatically invoked 
by REDUCE. For example, in the evaluation of the residuals of the 
nonlinear equations and in boundary conditions within the iteration. 

• £13-17 defines the operator psolv which is used to find the pressure 
correction through (|2^) . For each term in in the right-hand side, it 
integrates with the correct integration constant to give a contribution 
to the pressure. 
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£18-22 defines the operator usolv which is used to find the horizontal 
velocity correction through (^31). For each term of the form it 



a 






(1) 


m 




(9x™ 





solves an ode with the appropriate bottom and linearised free-surface 
boundary condition. 

£23-28 For compactness of the output, it is convenient to represent the 
film thickness and its derivatives via an operator h, £24. h is to depend 
upon X and t, £25, and h(ni) denotes d"^r]/dx"^, for example h(0) 
and h(l) denote rj and rj' respectively, and so x derivatives operate 
according to the transformations on £26-27. Whereas £28 encodes the 
fact that 

gmg 

where the dependence of rj upon x is recognised in the spatial deriva- 
tives of g. 

Note that there is no need to define a specific operator to compute the 
correction to the vertical velocity, from (^), because it is obtained simply 
by integrating the right-hand side as implemented by the standard int 
operator in REDUCE. 

The linear approximation is specified in £31. 

Iteration. The iteration is carried out until the errors are O {d"^ /dx"^), £35, 
whether in the form 77*™ or as 7]"ri^ or as some other such nonlinear combi- 
nation of derivatives of the film thickness rj. This is achieved by carrying a 
dummy variable d whose exponent counts the number of x derivatives, and 
discarding any term of order 7 or more in d. 

£36 establishes the expansion for r]"/ in a series in small slope r]'. 
Note the use of d to count the total number of x derivatives in each term. 

One variation in the iteration is that the loop, £37-55, uses new information 
as it becomes available: 

(a) first, the v correction from continuity, £38-40; 

(b) second, the pressure correction is found, £40-47, from the vertical mo- 
mentum residual and the normal stress across the free-surface (the 
term in v' on the left-hand side of is included in the right-hand 
side of £43 because we use the latest approximation to v); 

(c) third, the u correction is found, £48-53, from the horizontal momentum 
residual and the tangential stress across the free-surface; 

(d) and lastly, the latest version of the model evolution ^ = g from the 
kinematic free surface condition, £54. 
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Upon executing this program, I find not only the evolution equation (JT^), but 
also the velocity and pressure fields. Specifically, to errors of fifth-order in d/dx, 

u - {c-\e)v'v"\ 

V ~ V VV : 

^ , // 
p ^ —7] 

+lv'"v" - (1 + c) vvW" -{l+c- le) ri'r ■ 

Expressions such as these inform us of the dominant details of the fiow modelled 
by (p!^; higher-order terms were computed but I have not recorded them here. 



3.4 Aside 

Although this model describes the long-term dynamics of thin films, it is limited 
in its usefulness (even with gravitational effects included). For example, in the 
linearised problem ( pi]) at finite wavenumber, the leading branch of the spec- 
trum merges with the next branch, the gravest shear mode Ui = sin7rC/2, at a 
wavenumber k ^ 2 to form a pair of oscillatory decaying modes. Such oscilla- 
tions are the remnants, under the strong viscous dissipation in thin films, of the 
waves which surface tension can support. In many applications, such decaying 
waves seem important, for example see the review by ChangQ. However, the 
model ([T^) cannot describe the necessary oscillations because it has only one 
component and is only first-order in time. This appears at first sight to be a 
strong limitation to the practical usefulness of centre manifold techniques in ap- 
plications. But with some imagination we can modify the governing equations so 
that a centre manifold model is formed based on the two leading branches of the 
spectrum |Q. Such a model has much wider applicability because it is a much 



improved description at finite wavenumber and it resolves shorter transients in 
time. 



4 Conclusion 

In summary, the proposed algorithm for the computer algebra derivation of low- 
dimensional models of dynamical systems is relatively: 

• simple to implement, because the computation of the residual is via a direct 
coding of the governing differential equations; 

• reliable, because it relies upon the actual residual going to zero — any error 
is picked up by a lack of convergence; 
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• flexible, because it removes the explicit tyranny of primitive approaches, 
such as the method of multiple scales, in forcing highly specific scalings 
upon the parameters and variables in the model — instead centre manifold 
theory assures the model is accurate to the order of accuracy of the residual. 

Acknowledgements: I thank Valery Roy for stimulating discussions during 
the preparation of this work. 
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